package fem2.element.crack_sda2;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.math.LinAlg;

import java.util.ArrayList;

import math2.LinAlgUtilities;
import math2.MathUtilities;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.element.crack_sda.SeparationMode;
import fem2.enu.MeshPartType;
import fem2.enu.StressType;
import fem2.jlinalg.RotationMatrix;
import fem2.material.cohesive_law2.CohesiveLaw;

public class SOSConstantJumpTetrahedraElement extends ConstantJumpSDAElement {

	protected static final double sqrt_1_3 = Math.sqrt(1.0 / 3);

	protected SeparationMode bestMode = null;

	public SOSConstantJumpTetrahedraElement(MeshPart mp, MaterialModel mm, CohesiveLaw cl) {
		super(mp, mm, cl);
		setCrackMeshPartType(MeshPartType.TRI3);
		setIntegrationOrderOnCrack(1);
	}

	public SOSConstantJumpTetrahedraElement(MeshPart mp, MaterialModel mm, CohesiveLaw cl,
			StressType stressType) {
		this(mp, mm, cl);
		this.stressType = stressType;
	}

	@Override
	protected void computeEnhancedStrainOperator(MeshPart mp, IMatrix g, double... xi) {
		double coeff = -crackMeasure / elementMeasure;
		g.set(0, 0, coeff * normal[0]);
		g.set(1, 1, coeff * normal[1]);
		g.set(2, 2, coeff * normal[2]);

		g.set(3, 0, coeff * normal[1]);
		g.set(3, 1, coeff * normal[0]);

		g.set(4, 1, coeff * normal[2]);
		g.set(4, 2, coeff * normal[1]);

		g.set(5, 0, coeff * normal[2]);
		g.set(5, 2, coeff * normal[0]);

		IMatrix J = computeLocalProjectionOnCrack();
		IMatrix tmp = (IMatrix) g.clone();
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp, BLAM.NO_TRANSPOSE, J, 0.0, g);
	}

	@Override
	protected void computeProjectionOperator(MeshPart mp, IMatrix g, double... xi) {
		double coeff = -crackMeasure / elementMeasure;
		g.set(0, 0, coeff * normal[0]);
		g.set(1, 1, coeff * normal[1]);
		g.set(2, 2, coeff * normal[2]);

		g.set(3, 0, coeff * normal[1]);
		g.set(3, 1, coeff * normal[0]);

		g.set(4, 1, coeff * normal[2]);
		g.set(4, 2, coeff * normal[1]);

		g.set(5, 0, coeff * normal[2]);
		g.set(5, 2, coeff * normal[0]);
	}

	@Override
	protected boolean checkOnsetCriteria() {
		IVector sigma = getStress(0, stressType);
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		/*
		 * perform spectral decomposition
		 */
		IVector ps = new ArrayVector(dim);
		IMatrix pd = new Array2DMatrix(dim, dim);
		LinAlgUtilities.decomposePrincipalStress(dim, sigma, ps, pd);

		if (ps.get(2) < getCohesiveLaw().getTensileStrength()) {
			return false;
		} else {

			alpha = new ArrayVector(dim);

			/*
			 * compute B operator
			 */
			IMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
			double[] xi = getIntegrationPoint(0);
			computeB(mp, B, xi);

			/*
			 * check and filter all valid separation modes
			 */
			ArrayList<SeparationMode> validModes = new ArrayList<SeparationMode>();
			// int[][] configurations = new int[][] { { 0 }, { 1 }, { 2 }, { 3
			// }, { 0, 1 }, { 0, 2 },
			// { 0, 3 }, { 1, 2 }, { 1, 3 }, { 2, 3 } };
			int[][] configurations = new int[][] { { 0 }, { 1 }, { 2 }, { 3 }, { 0, 1 }, { 0, 2 },
					{ 0, 3 } };

			for (int i = 0; i < configurations.length; i++) {
				double gradx = 0, grady = 0, gradz = 0;
				for (int j = 0; j < configurations[i].length; j++) {
					gradx += B.get(0, configurations[i][j] * dim);
					grady += B.get(1, 1 + configurations[i][j] * dim);
					gradz += B.get(2, 2 + configurations[i][j] * dim);
				}
				double normgrad = Math.sqrt(gradx * gradx + grady * grady + gradz * gradz);
				gradx /= normgrad;
				grady /= normgrad;
				gradz /= normgrad;
				/*
				 * compute crack area
				 */
				double A = elementMeasure * normgrad;
				/*
				 * check the separation mode
				 */
				SeparationMode mode = new SeparationMode(mp, new ArrayVector(gradx, grady, gradz),
						A, configurations[i]);
				if (mode.validate()) {
					validModes.add(mode);
				}
			}

			/*
			 * realize crack normal
			 */
			double n1 = pd.get(0, 2);
			double n2 = pd.get(1, 2);
			double n3 = pd.get(2, 2);

			/*
			 * find all in valid separation modes, the mode make the smallest
			 * angle with crack normal vector from Rankine criterion
			 */
			double cos = 0.0;
			for (SeparationMode mode : validModes) {
				IVector N = mode.getSeparationNormal();
				double tmp = N.get(0) * n1 + N.get(1) * n2 + N.get(2) * n3;
				if (Math.abs(tmp) > Math.abs(cos)) {
					cos = tmp;
					bestMode = mode;
				}
			}

			/*
			 * get the normal from the best mode to become normal of the crack
			 */
			normal = ((ArrayVector) bestMode.getSeparationNormal()).getValues();
			crackMeasure = bestMode.getCrackArea();

			/*
			 * rotate the local crack frame obtain from Rankine criterion to the
			 * frame of best mode (strictly enforcing orthogonalization)
			 */
			double[] z = MathUtilities.cross(n1, n2, n3, normal[0], normal[1], normal[2]);
			double normz = LinAlg.norm2(dim, z);

			if (normz > MathUtilities.tol14) {
				/*
				 * find the angle between (n1,n2,n3) and (gradx,grady,gradz)
				 */
				cos = n1 * normal[0] + n2 * normal[1] + n3 * normal[2];
				double theta = Math.acos(cos);
				/*
				 * rotate other two principal stresses
				 */
				RotationMatrix rot = new RotationMatrix(0, 0, 0, z[0], z[1], z[2], theta);
				tangent1 = rot.timesXYZ(pd.get(0, 0), pd.get(1, 0), pd.get(2, 0));
				tangent2 = rot.timesXYZ(pd.get(0, 1), pd.get(1, 1), pd.get(2, 1));
			} else {
				tangent1 = new double[] { pd.get(0, 0), pd.get(1, 0), pd.get(2, 0) };
				tangent2 = new double[] { pd.get(0, 1), pd.get(1, 1), pd.get(2, 1) };
			}

			/*
			 * compute crack boundary
			 */
			crackBoundary = bestMode.getCrackBoundary();

			return true;
		}
	}

}
